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In a uniform fluid, a quantized vortex line with circulation h/M can support long-wavelength 
helical traveling waves oc e 1 ( kz ~ UJ k t ) w jth the well-known Kelvin dispersion relation ~ 
(hk 2 /2M) ln(l/|fc|£), where £ is the vortex-core radius. This result is extended to include the effect 
of a nonuniform harmonic trap potential, using a quantum generalization of the Biot-Savart law that 
determines the local velocity V of each element of the vortex line. The normal-mode eigenfunctions 
form an orthogonal Sturm-Liouville set. Although the line's curvature dominates the dynamics, the 
transverse and axial trapping potential also affect the normal modes of a straight vortex on the 
symmetry axis of an axisymmetric Thomas-Fermi condensate. The leading effect of the nonuniform 
condensate density is to increase the amplitude along the axis away from the trap center. Near the 
ends, however, a boundary layer forms to satisfy the natural Sturm-Liouville boundary conditions. 
For a given applied frequency, the next-order correction renormalizes the local wavenumber k(z) 
upward near the trap center, and k(z) then increases still more toward the ends. 

PACS numbers: 03.75.Kk, 67.40.Vs 



I. INTRODUCTION 

Quantized vortex lines constitute one of the most characteristic features of superfluidity. They were originally 
predicted to occur in superfluid He II 0, and have since been the subject of many studies [1- [j More recently, 
uantized vortices have been predicted and observed in the Bose-Einstein condensates of dilute trapped cold gases [a, 
21 HLIll' Apart from the quantized circulation k = h/M, their properties closely resemble those of classical vortex 
lines, which have long been known to support circularly polarized oscillatory normal modes |Toj . These Kelvin waves 
have a characteristic long- wavelength dispersion relation ujk ~ {nk^/Air) ln(l/fc£), where k is the wavenumber along 
the vortex axis, k is the circulation and £ is the vortex-core size. For a classical vortex with a hollow core, the general 
dispersion relation involves the Bessel functions Kq(Js£) and Kx(k£) 0,0, and the preceding approximate formula 
applies only for fc£ < 0.2. In this Kelvin mode, each element of the vortex core executes a circular orbit with a sense 
opposite to the circulating velocity around the vortex line. In the quantum-mechanical version, these modes therefore 
have quantized angular momentum with m — — 1 (in units of K) relative to the angular momentum of the vortex 
itself 020. 

In the last year, Bretin et al. have reported direct evidence for such Kelvin modes in an elongated cigar-shaped 
87 Rb condensate. By driving the condensate with a near-resonant quadrupole excitation (m = ±2), they observed 
preferential decay of the mode with m = — 2 relative to that for m = 2. The authors ascribed this difference to 
the creation of two Kelvin- wave excitations, each with m = — 1 and with opposite momenta ±hk determined by 
conservation of energy. Theoretical analysis confirms this process and the corresponding decay line shape |15|. In 
addition, transverse images of the expanded condensate show the periodic structure associated with the Kelvin-mode 
deformation [l4j . 

Essentially all theoretical analyses of Kelvin waves in a superfluid have assumed that the fluid is uniform along the 



naiysi 

vortex axis (see, however, Ref. [l6|) leading to the long- wavelength expression 

where the absolute value allows for the possibility of waves propagating in both directions. This unusual dispersion 
relation differs significantly from the nondispersive result u>k — ck for waves on a string. It arises because the dynamics 
here involves the velocity instead of the more familiar acceleration; in effect, a vortex is massless, and its motion obeys 
first-order dynamical equations 0, 0] . 
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In the case of a large dilute trapped Bose-Einstein condensate (BEC) in the Thomas- Fermi (TF) limit, the quantum 
of circulation h/M and the vortex core size £ = (87rart ) _1 ^ 2 are the only material parameters, where a is the s-wave 
scattering length and no is the central number density of the condensate. The simplest application is to a straight 
vortex initially at rest along the symmetry axis of a condensate. In this situation, it would seem natural to replace 
no by the ^-dependent TF density n(z) = no (l — z 2 / R 2 ) along the symmetry axis, where R z is the TF condensate 
dimension in the z direction. Such a substitution implies that the leading effect of the spatial inhomogeneity appears 
only in the logarithmic factor in Eq. Q . One aim of the present work is to demonstrate that there are other more 
significant alterations. 

A very general approach would be to rely on the time-dependent Gross-Pitaevskii (GP) equation 0, 0] that 
determines the time evolution of the condensate wave function ^ . Such a direct solution for a vortex wave in a 
nonuniform condensate would involve considerable numerical analysis. Another possibility is to generalize Pitaevskii's 
linearized theory of small-amplitude oscillations of a quantized vortex line |l2( to the case of an inhomogeneous 
condensate. This calculation rapidly becomes intricate, for it involves several very different length scales: £, fc _1 , R± 
and R z ; in contrast, Pitaevskii studied an unbounded condensate and required only the two parameters £ and fc _1 . 

A more promising procedure is to rely on the quantum analog of the Biot-Savart law that determines the local 
velocity of each element of the vortex |20j . This formula is readily specialized to the case of a nearly straight vortex 
line in a large TF condensate. The system can execute circularly polarized helical waves, and the associated normal- 
mode eigenfunctions constitute a Sturm-Liouville set subject to natural boundary conditions at the two ends. For 
most of the length, the spatial nonuniformity can be incorporated with a phase-integral (WKB) approximation, but 
a thin boundary layer forms near the ends. Section II provides the basic formalism, and the approximate solution 
is presented in Sec. III. Comparison with the experimental results and the limitations of the present analysis are 
discussed in Sec. IV. 



II. BASIC FORMALISM 



Consider an axisymmetric trap potential 



V tI (r,z) = -M(Lj 2 x r 2 +LJ 2 z z 2 ), 



(2) 



where r, <fi and z are cylindrical polar coordinates. In the TF limit of a large axisymmetric condensate, the particle 
density is given by 



(3) 



where R± = v / 2fi/Muj 2 J _ and R z = y/2fi/Mu) 2 are the radial and axial condensate dimensions, [i is the chemical 
potential and no = M/z/(47rft 2 a) is the central density. A detailed analysis of the time-dependent GP equation for a 
vortex in such a condensate yields the following expression [2(| for the velocity V(r, z) of an element of the vortex 
line located at (r, z): 



V(r,z) = 



2M 



-Ink' 





JC 2 \ 




Ik' 







i x W tr (r,z) 
fj,n(r,z)/n 



Kb 



(4) 



Here, t is the local tangent vector to the vortex axis, b is the binormal and K, is the curvature. For simplicity, the 
condensate is assumed to be stationary (nonrotating). 

In fact, the argument of the logarithm here requires some discussion. When Eq. is used to study the precession 
of a long straight vortex parallel to the symmetry axis of a disk-shaped condensate [21j, the curvature K. vanishes. 
Thus the logarithmic factor reduces to log(i?j_/£); it arises because the long-range circulating velocity field of the 
vortex line extends to infinity like 1 jr and is cut off only by the radius R± of the condensate. 

The situation is more complicated when the vortex undergoes a helical distortion like the Kelvin mode. The original 
derivation |lOj | of Eq. Q assumed an arbitrary initial axisymmetric azimuthal velocity Uo(r)4> for the unperturbed 
fluid. Fortunately, the situation becomes much simpler for the special case of a straight vortex line with a hollow 
core of radius £ in a uniform fluid. The unperturbed flow velocity Uo(r) — K/{2irr) is irrotational, and the resulting 
perturbed motion for the normal modes can now be characterized by a single scalar velocity potential <&(r,(j>,z,t) 
instead of the three separate components of the perturbations in the velocity. Indeed, the velocity potential for the 
Kelvin mode has the form 



<I>(r, (p, z, t) = $ exp [i (kz - cut - <j))\ Kx(kr), 



(5) 
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with angular momentum m = — 1. Here K\(kr) is a Bessel function that decays exponentially at infinity, as is 
appropriate for a solution of Laplace's equation that is periodic with wavenumber k along the z axis. In this picture |22| . 
the Kelvin mode is simply a surface wave propagating along the vortex core; the corresponding amplitude is confined to 
a cylindrical region with radius ~ k^ 1 . Consequently, for a propagating wave, the natural cutoff inside the logarithm 
in Eq. @ is the smaller of the two lengths R± and k~ l ] thus a useful approximation for a helical deformation is to 
replace the logarithm in Eq. J3J by 

C = ]n(^-=J==], (6) 

which holds with logarithmic accuracy (note that the curvature JC for a small-amplitude helical deformation is pro- 
portional to the radius of the deformation and to the squared wavenumber |23j). 

As additional justification for this approximation, it is helpful to recall the collective analog of Kelvin waves in 
a rotating vortex lattice with mean areal density n v = 2fl/n = MQ/ttH and intervortex separation b = {itn^ 1 / 2 . 
The relevant normal modes of this system are characterized by an axial wave number k and have a dispersion 
relation HUE HI 

2M \ki) 2 n ; {2n + (hk 2 /2M)\n(b/0, if fc& < 1. 

When kb is large, the wave amplitude arising from a given vortex vanishes exponentially at the position of the nearest 
neighbors. Thus the first and last terms in Eq. @ are negligible, and the dispersion relation reduces to that of a 
single vortex. In the opposite limit (kb <C 1), the leading (first) term w 2f2 is just the dispersion relation of a classical 
inertial wave in a uniformly rotating fluid [2jj , and the last term in Eq. (0 acts to alter the logarithmic factor. Once 
again the maximum length in the logarithm is the smaller of the relevant lengths (which here are b and k^ 1 ). 



A. Dynamical equations for a deformed vortex line 



The preceding equations can now be applied to an originally straight vortex along the axis of symmetry. The 
lateral displacement is specified by the two-dimensional vector x(z, £), y(z, t), and only terms of first order in x or y 
are retained. To this order, the tangent vector has the form 

i& z + x'x + y'y, (8) 

where the prime denotes a derivative with respect to z. Similarly, the gradient of the trap potential becomes 



(9) 



Finally, the curvature term becomes Kb w —y"x + x"y. A combination of these various terms with V(z, t) = xx + yy 
yields 



y 
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where C is defined in Eq. ©. 
Processing modes have the form 



l-z 2 /R 2 z \R 2 R\ 



x(z, t) = f(z) sinwt, 
y{z,t) = f(z)cosujt. 



(10) 

(11) 



(12) 
(13) 



Each element of the vortex core executes a clockwise circle (in the negative sense) with angular frequency u> (assuming 
uj > 0); for uj < 0, the sense of rotation reverses and becomes counterclockwise (in the positive sense). The unknown 
amplitude f(z) obeys the ordinary differential equation 



. nr. 

U)t = 

J 2M 



l-z 2 /R 2 z \ R 2 Z R 



(14) 
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Here, the term proportional to /' arises from the tipping of the vortex line and the trap confinement in the axial 
direction. The term proportional to / arises from the trap confinement in the radial direction; it occurs even for a 
straight vortex and accounts for the uniform precession of a straight vortex line in a trap. The term with /" arises 
from the curvature of the vortex line. 

In the limit of a purely two-dimensional radial trap with lu z — > and R z — ► oo, the solution is a uniform propagating 
wave / oc e lkz with the dispersion relation 

hk 2 i { R± 1 \ / 2 \ 

^ = w ln {- ^frm ) v - J ' (15) 

This expression is positive for k 2 R\ > 2, reproducing the Kelvin result QJ in the limit kR± 3> 1 because fc _1 then 
provides the relevant logarithmic cutoff. If, on the other hand, kR± is small, then the frequency becomes negative 
and motion approaches a uniform precession in the positive sense with angular frequency (h/MR\) ln(i?j_/£) |20|. as 
seen in recent experiments |28j . Note that this analysis predicts a stationary helical deformation when the curvature 
just balances the tendency to uniform precession (namely when k R\ = 2). 



B. Sturm-Liouville analysis of normal modes 



In the case of a nonuniform condensate with slowly varying density along the axis, the most interesting normal 
modes are those with many nodes, describable with a slowly varying wave number. In this case, the logarithmic 
parameter C in Eqs. © and (|14|l is approximately constant. It is convenient to introduce a dimensionless axial 
variable z with R z as length scale, so that — 1 < z < 1. 

The ordinary differential equation 114|) assumes the dimensionless Sturm-Liouville form 

Vf(z)-2^f(z) = T 2 (l-z 2 )f(z). (16) 



Here 



is a self-adjoint differential operator and 



r - hHl/ki) (18) 

is a large parameter. Note that T is approximately kR z , with k determined from the Kelvin dispersion relation for a 
uniform fluid (the experimental value k « 0.86 /im -1 evaluated with u> as half the quadrupole oscillation frequency 
for m = —2 gives T w 32). 

This equation closely resembles Legendre's equation [29l l3fl l3lj . and the physical solutions must obey similar 
"natural" boundary conditions at z = ±1 because of the factor (1 — z 2 ) in the operator T>. These conditions mean 
that f[z) must merely remain finite as z — ► ±1. A series expansion around z = leads to a three-term recursion 
relation, which makes such an approach quite intricate. In addition, the Frobenius method shows that z = ±1 arc 
regular singular points, with one finite solution in each case, as anticipated from the natural boundary conditions. 
The second solution near \z\ = 1 has a logarithmic singularity with leading term ~ ln(l — \z\). The simultaneous pair 
of finite- value conditions at both ends selects a discrete set of allowed values for T 2 , which is effectively the eigenvalue. 

Let j be the discrete index for the allowed solutions. The real eigenfunctions fj(z) obey the eigenvalue problem 

Vfi{z) - 2§k- W = T 2 (1 - z 2 ) f 3 {z) (19) 
on the interval — 1 < z < 1, with the eigenvalue V 2 . The solutions are orthogonal with 

dzMz)(l-z 2 )f 3 (z) = Xijtj. (20) 
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III. CONSTRUCTION OF APPROXIMATE EIGENFUNCTIONS AND EIGENVALUES 

To construct the approximate solutions of Eq. (|19fl . it is necessary to approach the problem in two separate but 
overlapping regions. Away from the ends, the WKB method gives a simple description, but the resulting "outer" 
solutions do not satisfy the natural boundary condition. Instead, the exact solution develops a thin boundary layer (the 
"inner" solution) that requires a different method. Matching the two solutions in the region of common applicability 
provides the required eigenvalue condition |32j| . 

A. Phase-integral (WKB) approximation 

The basic problem is to find an approximate solution of Eq. (|19|l that explicitly includes the density variation that 
leads to the factors 1 — z 2 . For an off-center straight vortex in a disk-shaped condensate, the radial analog of this 
factor leads to a precession frequency that increases with increasing vortex distance from the axis of symmetry 0, l33| . 
Here, in contrast, the density varies along the direction of the helical wave. The essential idea is that the properties of 
the medium vary only on the overall length scale R z . Except for the lowest normal modes, one can therefore define a 
local axial wavenumber k, with kR z 3> 1. Based on the experiment of Bretin et al. |l4j . it is convenient to impose the 
additional restriction kR± 3> 1 (the experimental value is kR± w 4), and to assume that fc£ is small [the experimental 
value is k£ w |, so that ln(l/fc£) w 1.4 is not especially large]. 

The phase- integral (WKB) method provides a natural approach to this problem [3^ . l34l l35| . This method assumes 
that 

f(z) = A(z)exp[irS(z)}. (21) 
Similar to the familiar quantum- mechanical case [3(| , the real and imaginary parts of Eq. I|19(l yield 

T 2 A = A' - — |f-A + T 2 S' 2 A ~ A", (22) 

1 — z 2 1 — z z 

= -^%S'A-2S'A'-S"A. (23) 
1 — z z 

The dominant term arises from the phase S, and it is useful to assume the following expansion in inverse powers of T: 

S = S + T- 2 Si + • • • , (24) 
A = A +T- 2 A 1 + --- . (25) 



To leading order in T, the real part (|22|l gives 



so that 



S * = 1, (26) 



S (z) = ±z. (27) 



In this order, the phase is simply that of a Kelvin wave [in dimensional units, So — kz, where k is determined as the 
solution of Eq. JTJ]. 

The next-order contribution is from the imaginary part. Since Sq = 0, Eq. immediately yields 

;A Q = A' , (28) 



1 

which reflects the balance between the longitudinal trap confinement on the left and the curvature on the right. Direct 
integration gives 

Ao(z) = -j=L=. (29) 

To this approximation, the solution resembles the usual WKB solution j3|| 

1 [ cos(rz), for even solutions; 
u (z) w { ( 30 ) 

VI — z 2 I sin(rz), for odd solutions. 

The leading effect of the axial variation caused by the nonuniform trap potential is to increase the amplitude of the 
normal mode as \z\ increases away from the origin, while the phase increases uniformly with z. 



6 



B. Structure of boundary layer near z = ±1 

The outer (WKB) solutions l|30[l apparently diverge for \z\ — > 1, which violates the natural boundary condition 
that the amplitude remain finite. Thus the exact solution must change character within a thin boundary layer of 
dimensionless thickness 5 <C 1 to be determined below. For definiteness, consider the positive end z = 1. It is 
convenient to introduce an "inner" variable Z by the relation 

Z =^T> (31) 

which gives an appropriate boundary layer near the surface at z = 1. 

Rewrite the eigenfunction f(z) as F(Z) in terms of the new variable, using z = 1 — ZS, and substitute into Eq. H16[l . 
This process yields an ordinary differential equation 

(2Z - Z 2 5) ^-J + 2 (1 - ZS) % + 2^F + T 2 S 2 (2Z - Z 2 5) F = 0. (32) 
cLZj " clZj 

Assume an expansion 

F{Z) ^ F (Z) + 5F 1 (Z) + ■ ■ ■ , (33) 
and substitute into Eq. ll-'il't . For small 5, the leading contribution is 



d 2 F 1 dF 



rVF = 0. (34) 



dZ 2 Z dZ 

It is convenient to make a specific choice of scale factor and take 

5 = i < 1 (35) 

as the dimensionless boundary-layer thickness. In this case, Eq. (|34|l reduces to Bessel's equation of order zero, with 
linearly independent solutions Jo(Z) and Y (Z). The latter is logarithmically singular as Z — * and must be rejected, 
leaving the inner solution that obeys the natural boundary condition 

F (Z) = C J (Z), (36) 

where Cq is a constant. 

Away from the surface (Z> 1), the boundary-layer solution 1)36(1 has the asymptotic form 

F {Z)~cJ^- (37) 



ttZ V 4 



For comparison, the WKB solution l|3U|) can also be expressed in terms of the inner variable Z, using y/l-z 2 w \/2Z5 
for <5 — > 0. For the even solution, this procedure yields 

m - cos( ^— z) . (38) 

Comparison of the overall factors gives the coefficient 



Co = ±^f (39) 
using T6 = 1. Matching the phase of the two oscillatory solutions fixes the eigenvalue 

Tj w M + ^ J 7r for the even solutions, (40) 

where Co is positive (negative) for j even (odd) . Similarly, the eigenvalue for the odd solutions is 

Tj » ( j + ^ J 7r for the odd solutions. (41) 
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FIG. 1: Typical matched solution for the eigenfunctions corresponding to large excitation quantum number (here for j = 10). 
The solid line is the outer WKB solution valid away from the ends, and the dashed line is the inner boundary-layer solution. 
Curve (a) is the even solution with F m 10.257T and curve (b) is the corresponding odd solution with V ~ 10.757T (both shown 
only for z > 0). 



Figure 1 shows the even and odd solutions for j = 10. The solid lines are the WKB approximations and the dashed 
lines are the boundary-layer solutions. Apart from the small discrepancy near the local maxima and minima, the fit 
is excellent. In practice, it suffices to take the inner solution Fq(Z) for the interval < Z < a\ where a\ « 2.4048 
is the first zero of Jo{Z). For the remaining interval < z < 1 — Sax, the leading WKB solution fo(z) gives a good 
description. 

The normalization integral 

1 = J ^dzf j (z) 2 (l-z 2 )=2^ dzf 3 (zf(l~z 2 ) (42) 
can be evaluated using the same approximations. Write I = I\ + 1%, where 

pza pi 

h = 2 dz fj(z) 2 (l — z 2 ) and I 2 = 2 dz fj(z) 2 (l - z 2 ) , (43) 

JO Jzo 

with zq = 1 — Sa\. 

For Ix, use of the WKB approximation immediately gives I\ w 1 as 6 — > 0. For I2, the boundary-layer solution 
F (Z) readily leads to 

I 2 ^25 f 1 dZ [Co Jo (Z)} 2 (2ZS) , (44) 
Jo 

which vanishes in the limit 5 — > (note that Cq oc <5 _1 ). 

The experiment [T3 | effectively drives the vortex with a quadrupole perturbation of given frequency. Depending on 
the width of the various resonances, this driving force will excite one or more of the normal modes, and the experiment 
seems to observe a nearly pure mode. The appropriate j value would be determined from conservation of energy and 
the symmetry of the normal modes. Unfortunately, a more direct observational study of the form of the eigenfunctions 
seems impractical. 
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C. Next-order WKB description 



It is interesting to consider the next-order contribution for large T, which arises from the first correction to the real 
part iPUl: 



R 2 



Q / ni _ A'q _ 2z A' 

~ A 1 - z 2 Ao ^ R2 I 



(45) 



where the last term reflects both the transverse and the axial trap confinement. For a wave with S' = 1, the explicit 
form of Aq leads to 



S[(z) 



R~ 



1 



1 



R\ 1 - z 2 



2(1 -z 2 ) 2 ' 

and both terms act to increase the phase because S[ is positive. Direct integration gives 



(46) 



R\ 1 . , 
i# + 4 lln 



1 + z 



1 



21- z 2 



(47) 



A combination with Eq. i|30|) now includes the effect of the nonuniformity on both the phase and the amplitude 



/(*) 



1 



cos 
sin 



Tz + ±- 
2T 



^ + 4 lln 



1 + z 
1 - z 



1 



21- z 2 



An expansion of the phase for small \z\ immediately yields the effective dimensionless wave number 

r ^ r+ r(lM)> r < 



(48) 



(49) 



where the first term in the correction arises from the overall trap confinement and the second term arises from the 
combined effect of the curvature and the longitudinal confinement. Note that r e ff exceeds T obtained as the solution 
of Eq. 0. 

More generally, the quantity S'(z) defines the effective dimensional wavenumber k e s(z) through the relation 



S'(z) = k cS (z)R z « kR z 



1 



1 



1 



1 



f 



k 2 R\ 1 — z 2 2k 2 R 2 (I - z 2 ) 2 



(50) 



where k is the solution of Eq. (ignoring the slow variation of the logarithmic factor). For small \z\, this expression 
reduces to (|49|l . but it now includes the region where I — z 2 is not small. Eventually this approach will fail near 
the ends of the condensate, but it does indicate that the Kelvin wave in a TF condensate should have not only an 
increased amplitude but also an increased wavenumber as \z\ increases from the center of the trap along the symmetry 
axis. 



IV. DISCUSSION AND POTENTIAL COMPARISON WITH EXPERIMENT 

As noted in Ref. the Kelvin dispersion relation Q for the ENS experiments predicts a wavenumber k ss 
0.86 /zm . This value agrees well with the observations after accounting for the expansion of the condensate (which 
is strong in the radial direction but only weak in the axial direction for the ENS trap, where R z /R± ~ 8.3). Since 
r w 32, the corrections studied in the present paper are formally small (of order T~ 2 ). Note, however, the important 
effect of trap anisotropy through the factor R 2 z jR\ ~ 70 in Eq. (|49|l . As a result, the correction is instead of order 
(kR±)~ 2 ~ 0.07, which might well be observable. 

Unfortunately, any direct comparison with Ref. [T^j is suspect because the original unperturbed vortex was bent. 
Such behavior was predicted in |38t 13^. l40| and subsequently observed in In contrast, the present theoretical 

description assumes a straight unperturbed vortex line along the axis of symmetry, and an extension to an initially 
bent vortex line would be difficult. 

A straight vortex line has been proved to be the ground state for a spherical and a flattened disk-shaped trap with 
Rz/R± < 1 0- More detailed numerical work I4p fi nds that the bent vortex first becomes the stable state for 
Rz/R^l ^ 2, which agrees with the observations in [Uj where R z /R± ~ 8.2. Thus it would be interesting to study 
Kelvin waves in a condensate that is only slightly elongated, in which case the present theoretical description should 
apply. 



9 

Acknowledgments 

I am grateful to J. Dalibard, S. Stringari, and A. Svidzinsky for stimulating discussions that initiated this work. 



L. Onsager, Nuovo Cimento 6, Suppl. 2, 249 (1949). 

R. P. Feynman, in Progress in Low Temperature Physics, edited by C. J. Gorter (North-Holland, Amsterdam, 1955), 
Vol. 1, Chap. 2. 

A. L. Fetter, in The Physics of Liquid and Solid Helium, edited by K. H. Bennemann and J. B. Ketterson (John Wiley 
and Sons, New York, 1976), Part I, pp. 207-305. 

R. J. Donnelly, Quantized Vortices in Helium II (Cambridge University Press, Cambridge, 1991), Chap. 2. 

M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. Hall, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 83, 2498 

(1999). 

K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Phys. Rev. Lett. 84, 806 (2000). 

J. R. Abo-Shaeer, C. Raman, J. M. Vogels, and W. Ketterle, Science 292, 476 (2001). 

P. C. Haljan, I. Coddington, P. Engels, and E. A. Cornell, Phys. Rev. Lett. 87, 210403 (2001). 

A. L. Fetter and A. A. Svidzinsky, J. Phys.: Condens. Matter 13, R135 (2001). 

W. Thomson (Lord Kelvin), Philos. Mag. 10, 155 (1880). 

Reference Q, pp. 79 and 179-181. 

L. P. Pitaevskii, Sov. Phys. JETP 13, 451 (1961). 

T. Isoshima and K. Machida, Phys. Rev. A 59, 2203 (1999). 

V. Brctin, P. Rosenbusch, F. Chevy, G. V. Shlyapnikov, and J. Dalibard, Phys. Rev. Lett. 90, 100403 (2003). 
T. Mizushima, M. Ichioka, and K. Machida, Phys. Rev. Lett. 90, 180401 (2003). 

F. Chevy and S. Stringari, Phys. Rev. A 68, 053601 (2003), take a different approach to study the collective modes of a 
dense vortex array in a trap, using diffused vorticity and rotational hydrodynamics to study the normal modes. 
A. L. Fetter and P. C. Hohenberg, Phys. Rev. 159, 330 (1967). 

D. Stauffer and A. L. Fetter, Phys. Rev. 168, 156 (1968). 

E. P. Gross, Nuovo Cimento 20, 454 (1961). 

A. A. Svidzinsky and A. L. Fetter, Phys. Rev. A 62, 063617 (2000). 
A. A. Svidzinsky and A. L. Fetter, Phys. Rev. Lett. 84, 5919 (2000). 

A. L. Fetter, in Lectures in Theoretical Physics, edited by K. Mahanthappa and W. E. Brittin (Gordon and Breach, N. Y., 
1969), Vol. X1B, pp. 337-344. 

J. S. R. Chisholm and R. M. Morris, Mathematical Methods in Physics (Saunders, Philadelphia, 1965), Sec. 11.3. 
E. S. Raja Gopal, Ann. Phys. (NY.) 29, 350 (1964). 

A. L. Fetter^ Phys. Rev. 162, 143 (1967). 
Reference Q, Sec. 6.2. 

S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability (Dover, New York, 1981), Sees. 23 and 68. 

B. P. Anderson, P. C. Haljan, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 88, 2857 (2000). 

R. Courant and D. Hilbert, Methods of Mathematical Physics (Interscience, New York, 1953), Vol. 1, pp. 208-211 and 
325-327. 

N. N. Lebedev, Special Functions and Their Applications (Dover, New York, 1972), Chaps. 4 and 7. 

A. L. Fetter and J. D. Walecka, Theoretical Mechanics of Particles and Continua (Dover, New York, 2003), Sec. 40 and 
App. D2. 

C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers (McGraw-Hill, New York, 
1978), Chaps. 9 and 10. 

E. Lundh and P. Ao, Phys. Rev. A 61, 063612 (2000). 

J. Heading, An Introduction to Phase-Integral Methods (Methuen, London, 1962), Chap. 2. 

The WKB approximation has also been used to treat the Bogoliubov equations in A. Csordas, R. Graham, and P. Szepfalusy, 
Phys. Rev. A 57, 4669 (1998). 

See, for ex amp le, D. J. Griffiths, Introduction to Quantum Mechanics (Prentice-Hall, Englewood Cliffs, NJ, 1995), Chap. 8. 
Reference |30|. Chaps. 5 and 6. 

J. J. Garci'a-Ripoll and V. M. Perez-Garda, Phys. Rev. A 63, 041603(R) (2001). 
J. J. Garci'a-Ripoll and V. M. Perez- Garcia, Phys. Rev. A 64, 053611 (2001). 
A. Aftalion and T. Riviere, Phys. Rev. A 64, 043611 (2001). 
P. Rosenbusch, V. Bretin, and J. Dalibard, Phys. Rev. Lett. 89, 200403 (2002). 
A. Aftalion and R. L. Jerrard, Phys. Rev. A 66, 023611 (2002). 

M. Modugno, L. Pricoupenko, and Y. Castin, Eur. Phys. J. D 22, 235 (2003), see Fig. 7. 



